      integer function getlin(line,ch)
      character line(*)
      integer ch,i,j,fgetc
      i=1
23000 continue
      j=fgetc(ch,line(i))
      if(j .ne. 0)then
      line(i)=char(0)
      getlin=j
      return
      endif
      if(line(i) .eq. char(10))then
      line(i)=char(0)
      getlin=0
      return
      endif
      i=i+1
23001 goto 23000
23002 continue
      end
      subroutine putlin(line,ch)
      integer ch,i,j,fputc
      character line(*)
      i=1
23007 if(.not.(line(i) .ne. char(0)))goto 23009
      j=fputc(ch,line(i))
23008 i=i+1
      goto 23007
23009 continue
      return
      end
      subroutine error(str)
      character str(*)
      call putc(char(10))
      call putlin(str,6)
      stop
      end
      subroutine putdec(cha,n,w)
      character ch(79)
      integer itoc,i,j,n,nd,w,cha,fputc
      nd=itoc(n,ch,79)
      i=nd+1
23010 if(.not.(i.le.w))goto 23012
      j=fputc(cha,' ')
23011 i=i+1
      goto 23010
23012 continue
      i=1
23013 if(.not.(i.le.nd))goto 23015
      j=fputc(cha,ch(i))
23014 i=i+1
      goto 23013
23015 continue
      return
      end
      integer function itoc(int,str,size)
      integer iabs, mod, d,i,int,intval,j,size
      character k, str(size), digits(10)
      data digits /'0','1','2','3','4','5','6','7','8','9'/
      intval=iabs(int)
      str(1)=char(0)
      i=1
23016 continue
      i=i+1
      d=mod(intval,10)
      str(i)=digits(d+1)
      intval=intval/10
23017 if(.not.(intval.eq.0.or.i.ge.size))goto 23016
23018 continue
      if(int.lt.0 .and. i.lt.size)then
      i=i+1
      str(i)='-'
      endif
      itoc=i-1
      j=1
23021 if(.not.(j.lt.i))goto 23023
      k=str(i)
      str(i)=str(j)
      str(j)=k
      i=i-1
23022 j=j+1
      goto 23021
23023 continue
      return
      end
      integer function ctoi(in,i)
      character in(1)
      character*10 digits
      integer index,d,i
      data digits/'0123456789'/
23024 if(in(i).eq.' ')then
      i=i+1
      goto 23024
      endif
23025 continue
      ctoi=0
23026 if(.not.(in(i).ne.char(0)))goto 23028
      d=index(digits,in(i))
      if(d.eq.0)then
      goto 23028
      endif
      ctoi=10*ctoi + d - 1
23027 i=i+1
      goto 23026
23028 continue
      return
      end
      subroutine mpcopy(a,b)
      integer a(1), b(1), i, d
      d=a(1)
      do23031 i=1,d
      b(i)=a(i)
23031 continue
23032 continue
      return
      end
      subroutine add(a,b,c)
      integer a(1), b(1), c(1), carry, i, t, ba
      data ba/1073741824/
      if(a(1).eq.b(1))then
      carry=0
      i=2
23035 if(.not.(i.le.a(1)))goto 23037
      t = a(i) + b(i) + carry
      if(t.lt.ba)then
      carry=0
      c(i)=t
      else
      carry=1
      c(i)=t-ba
      endif
23036 i=i+1
      goto 23035
23037 continue
      if(carry.eq.1)then
      c(1)=a(1)+1
      i=c(1)
      c(i)=1
      else
      c(1)=a(1)
      endif
      else
      if(a(1).lt.b(1))then
      call addl(a,b,c)
      else
      call addl(b,a,c)
      endif
      endif
      return
      end
      subroutine addl(a,b,c)
      integer a(1), b(1), c(1), carry, i, t, ba
      data ba/1073741824/
      carry=0
      i=2
23044 if(.not.(i.le.a(1)))goto 23046
      t = a(i) + b(i) + carry
      if(t.lt.ba)then
      carry=0
      c(i)=t
      else
      carry=1
      c(i)=t-ba
      endif
23045 i=i+1
      goto 23044
23046 continue
23049 if(.not.(i.le.b(1)))goto 23051
      t = b(i) + carry
      if(t.lt.ba)then
      carry=0
      c(i)=t
      else
      carry=1
      c(i)=t-ba
      endif
23050 i=i+1
      goto 23049
23051 continue
      if(carry.eq.1)then
      c(1)=b(1)+1
      i=c(1)
      c(i)=1
      else
      c(1)=b(1)
      endif
      return
      end
      subroutine sub(a,b,c)
      integer a(1), b(1), c(1), borrow, i, t, ba
      character*80 smsg
      data ba/1073741824/
      smsg = 'negative diff in sub'//char(0)
      if(a(1) .lt. b(1))then
      call error(smsg)
      endif
      borrow=0
      i=2
23058 if(.not.(i.le.b(1)))goto 23060
      t = a(i) - b(i) - borrow
      if(t.lt.0)then
      borrow=1
      c(i)=t+ba
      else
      borrow=0
      c(i)=t
      endif
23059 i=i+1
      goto 23058
23060 continue
23063 if(.not.(i.le.a(1)))goto 23065
      t = a(i) - borrow
      if(t.lt.0)then
      borrow=1
      c(i)=t+ba
      else
      borrow=0
      c(i)=t
      endif
23064 i=i+1
      goto 23063
23065 continue
      if(borrow.eq.1)then
      call error(smsg)
      else
      c(1)=a(1)
      i=c(1)
23070 if(i.gt.2 .and. c(i).eq.0)then
      c(1)=c(1)-1
      i=c(1)
      goto 23070
      endif
23071 continue
      endif
      return
      end
      subroutine mul1(a,b,c)
      integer a(1), b, c(1), carry, d, i, x(2)
      character*80 negb
      negb = ' negative b in mul1'//char(0)
      if(b.gt.0)then
      goto 10
      endif
      if(b.lt.0)then
      call putdec(6,401,6)
      call putdec(6,b,10)
      call putlin(negb,6)
      call putc(char(10))
      endif
      c(1)=2
      c(2)=0
      return
10    carry=0
      d=a(1)
      do23076 i=2,d 
      call muladd(a(i),b,carry,x)
      carry = x(2)
      c(i) = x(1)
23076 continue
23077 continue
      if(carry.eq.0)then
      c(1)=d
      else
      c(1)=d+1
      i=c(1)
      c(i)=carry
      endif
      return
      end
      integer function div1(a,b,q)
      integer a(1), b, q(1), i, r, y(2)
      character*80 negb
      negb = '  b <= 0 in div1'//char(0)
      if(b.gt.0)then
      goto 10
      endif
      call putdec(6,501,6)
      call putdec(6,b,10)
      call putlin(negb,6)
      call putc(char(10))
      stop
10    r=0
      i=a(1)
23082 if(.not.(i.gt.1))goto 23084
      call divrem(r,a(i),b,y)
      r=y(2)
      q(i)=y(1)
23083 i=i-1
      goto 23082
23084 continue
      div1=r
      q(1)=a(1)
      i=q(1)
23085 if(i.gt.2 .and. q(i).eq.0)then
      q(1)=q(1)-1
      i=q(1)
      goto 23085
      endif
23086 continue
      return
      end
      subroutine mpchk(code,a)
      integer code, a(1), d, i
      character*80 s1,s2,s3,s4,s5
      s5 = ' mpchk'//char(0)
      s1 = ' dimension < 2 '//char(0)
      s2 = ' dimension > MPDIM '//char(0)
      s3 = ' negative digit '//char(0)
      s4 = ' digit >= MPBASE '//char(0)
      d=a(1)
      call putlin(s5,6)
      call putdec(6,code,5)
      do23087 i=1,d
      call putdec(6,a(i),10)
23087 continue
23088 continue
      call putc(char(10))
      if(d.lt.2)then
      goto 10
      endif
      if(d.gt.120)then
      goto 20
      endif
      do23093 i=2,d 
      if(a(i).lt.0)then
      goto 30
      endif
      if(a(i).ge.1073741824)then
      goto 40
      endif
23093 continue
23094 continue
      return
10    call putdec(6,code,5)
      call putdec(6,d,10)
      call putlin(s1,6)
      goto 90
20    call putdec(6,code,5)
      call putdec(6,d,10)
      call putlin(s2,6)
      goto 90
30    call putdec(6,code,5)
      call putdec(6,d,10)
      call putdec(6,i,10)
      call putdec(6,a(i),10)
      call putlin(s3,6)
      goto 90
40    call putdec(6,code,5)
      call putdec(6,d,10)
      call putdec(6,i,10)
      call putdec(6,a(i),10)
      call putlin(s4,6)
90    call putc(char(10))
      stop
      end
      subroutine mul(a,b,c)
      integer a(1), b(1), c(1), i,da,db,j,carry,x(2),k
      if(a(1).lt.3)then
      call mul1(b,a(2),c)
      return
      endif
      if(b(1).lt.3)then
      call mul1(a,b(2),c)
      return
      endif
      da=a(1)
      db=b(1)
      do23103 i=2,da
      c(i)=0
23103 continue
23104 continue
      do23105 j=2,db 
      carry=0
      do23107 i=2,da 
      k=i+j-2
      call muladd(a(i),b(j),c(k)+carry,x)
      carry=x(2)
      c(k)=x(1)
23107 continue
23108 continue
      k=da+j-1
      c(k)=carry
23105 continue
23106 continue
      c(1)=da+db-1
      i=c(1)
      if(c(i).eq.0)then
      c(1) = c(1) - 1
      endif
      return
      end
      subroutine div(a,b,q,r)
      integer a(1),b(1),q(1),r(1),m,d,db,j,qhat,l,i
      integer div1,y(2),x(2),borrow,tmp,k
      common /mpcomm/ junkc(600),junk1(120),junk2(120),junk3(120),junk4(
     x120)
      character junkc
      integer junk1,junk2,junk3,junk4
      character*80 norm,q1,unn
      norm = 'div norm failure'//char(0)
      q1 = ' qhat - 1 '//char(0)
      unn = 'div unnorm failure'//char(0)
      if(b(1).lt.3)then
      r(1)=2
      r(2)=div1(a,b(2),q)
      return
      endif
      m=a(1)-b(1)+2
      db=b(1)
      if(m.gt.1)then
      d=1073741824/(1+b(db))
      call mul1(a,d,junk1)
      call mul1(b,d,junk2)
      k=junk2(1)
      junk2(k+1)=0
      if(k.ne.db)then
      call putdec(6,db,10)
      call putdec(6,k,10)
      call putc(char(10))
      call error(norm)
      endif
      l=junk1(1)
      if(l.eq.a(1))then
      l=l+1
      junk1(l)=0
      endif
      j=0
23119 if(.not.(j.lt.m-1))goto 23121
      if(junk1(l).eq.junk2(db))then
      qhat=1073741824-1
      i=junk1(l-1)+junk2(db)
      else
      call divrem(junk1(l),junk1(l-1),junk2(db),y)
      i=y(2)
      qhat=y(1)
      endif
      call muladd(qhat,junk2(db-1),0,x)
23124 if(x(2).gt.i.or.(x(2).eq.i.and.x(1).gt.junk1(l-2)))then
      qhat=qhat-1
      i=i+junk2(db)
      call muladd(qhat,junk2(db-1),0,x)
      goto 23124
      endif
23125 continue
      call mul1(junk2,qhat,junk3)
      if(junk2(1).eq.junk3(1))then
      k=junk2(1)+1
      junk3(1)=k
      junk3(k)=0
      endif
      borrow=0
      k=l-junk2(1)+1
      i=2
23128 if(.not.(i.le.junk3(1)))goto 23130
      tmp=junk1(k)-junk3(i)-borrow
      if(tmp.lt.0)then
      borrow=1
      junk1(k)=tmp+1073741824
      else
      borrow=0
      junk1(k)=tmp
      endif
      k=k+1
23129 i=i+1
      goto 23128
23130 continue
      k=m-j
      if(borrow.eq.0)then
      q(k)=qhat
      else
      q(k)=qhat-1
      borrow=0
      k=l-junk2(1)+1
      i=2
23135 if(.not.(i.le.junk2(1)))goto 23137
      tmp=junk1(k)+junk2(i)+borrow
      if(tmp.lt.1073741824)then
      borrow=0
      junk1(k)=tmp
      else
      borrow=1
      junk1(k)=tmp-1073741824
      endif
      k=k+1
23136 i=i+1
      goto 23135
23137 continue
      endif
      l=l-1
23120 j=j+1
      goto 23119
23121 continue
23140 if(db.gt.2.and.junk1(db).eq.0)then
      db=db-1
      goto 23140
      endif
23141 continue
      junk1(1)=db
      i=div1(junk1,d,r)
      if(i.ne.0)then
      call putdec(6,i,10)
      call putc(char(10))
      call error(unn)
      endif
      if(q(m).eq.0 .and. m.gt.2)then
      q(1)=m-1
      else
      q(1)=m
      endif
      else
      call mpcopy(a,r)
      q(1)=2
      q(2)=0
      endif
      return
      end
      integer function mpitoc(int,str,size)
      integer int(1),size,div1,d,i,j
      character str(1),digits(10),k
      character*80 big
      common /mpcomm/ junkc(600),junk1(120),junk2(120),junk3(120),junk4(
     x120)
      character junkc
      integer junk1,junk2,junk3,junk4
      data digits/'0','1','2','3','4','5','6','7','8','9'/
      big = 'int too big for mpitoc'//char(0)
      if(int(1) .gt. 120)then
      call error(big)
      endif
      call mpcopy(int,junk1)
      str(1)=char(0)
      i=1
23148 continue
      i=i+1
      d=div1(junk1,10,junk1)
      str(i)=digits(d+1)
23149 if(.not.((junk1(1).eq.2.and.junk1(2).eq.0).or.i.ge.size))goto 2314
     x8
23150 continue
      mpitoc=i-1
      j=1
23151 if(.not.(j.lt.i))goto 23153
      k=str(i)
      str(i)=str(j)
      str(j)=k
      i=i-1
23152 j=j+1
      goto 23151
23153 continue
      return
      end
      subroutine mpctoi(in,i,out)
      character in(1)
      integer index, d, i, out(1), x(2)
      equivalence (x(2),d)
      character*10 digits
      data digits /'0123456789'/
      data x(1)/2/
23154 if(in(i).eq.' ')then
      i=i+1
      goto 23154
      endif
23155 continue
      out(1)=2
      out(2)=0
23156 if(.not.(in(i).ne.char(0)))goto 23158
      d=index(digits,in(i))
      if(d.eq.0)then
      goto 23158
      endif
      call mul1(out,10,out)
      d=d-1
      call add(out,x,out)
23157 i=i+1
      goto 23156
23158 continue
      return
      end
      subroutine mptdec(lunit,n,w)
      integer mpitoc,i,n(1),nd,w,lunit,j,fputc
      common /mpcomm/ junkc(600),junk1(120),junk2(120),junk3(120),junk4(
     x120)
      character junkc
      integer junk1,junk2,junk3,junk4
      nd=mpitoc(n,junkc,600)
      i=nd+1
23161 if(.not.(i.le.w))goto 23163
      j= fputc(lunit,' ')
23162 i=i+1
      goto 23161
23163 continue
      i=1
23164 if(.not.(i.le.nd))goto 23166
      j= fputc(lunit,junkc(i))
23165 i=i+1
      goto 23164
23166 continue
      return
      end
      integer function mpcmp(a,b)
      integer a(1), b(1), i
      if(a(1).lt.b(1))then
      mpcmp=-1
      else
      if(a(1).gt.b(1))then
      mpcmp=1
      else
      mpcmp=0
      i=a(1)
23171 if(.not.(i.gt.1))goto 23173
      if(a(i).lt.b(i))then
      mpcmp=-1
      goto 23173
      else
      if(a(i).gt.b(i))then
      mpcmp=1
      goto 23173
      endif
      endif
23172 i=i-1
      goto 23171
23173 continue
      endif
      endif
      return
      end
      integer function mpodiv(a)
      integer a(1),div1,i,d,j,l
      character*80 e
      e = 'error in MPODIV'//char(0)
      if(a(1).eq.2.and.a(2).eq.0)then
      mpodiv=-1
      return
      endif
      i=2
23180 if(.not.(a(i).eq.0))goto 23182
23181 i=i+1
      goto 23180
23182 continue
      i=i-2
      if(i.gt.0)then
      d=a(1)-i
      do23185 j=2,d 
      l=j+i
      a(j)=a(l)
23185 continue
23186 continue
      a(1)=d
      endif
      j=0
      d=a(2)
23187 if(mod(d,2).eq.0)then
      j=j+1
      d=d/2
      goto 23187
      endif
23188 continue
      mpodiv=30*i+j
      if(j.gt.0)then
      if(div1(a,2**j,a).ne.0)then
      call mpchk(550,a)
      call error(e)
      endif
      endif
      return
      end
      subroutine mpgocd(a,b,c)
      integer a(1),b(1),c(1),i,j,mpodiv
      common /mpcomm/ junkc(600),junk1(120),junk2(120),junk3(120),junk4(
     x120)
      character junkc
      integer junk1,junk2,junk3,junk4
      call mpcopy(a,c)
      call mpcopy(b,junk1)
      i=mpodiv(c)
      if(mpodiv(junk1).lt.0)then
      return
      endif
      if(i.lt.0)then
      call mpcopy(junk1,c)
      return
      endif
      i=mpcmp(c,junk1)
23197 if(i.ne.0)then
      if(i.lt.0)then
      call sub(junk1,c,junk1)
      j=mpodiv(junk1)
      else
      call sub(c,junk1,c)
      j=mpodiv(c)
      endif
      i=mpcmp(c,junk1)
      goto 23197
      endif
23198 continue
      return
      end
      subroutine mpower(b,e,m,a)
      integer b(1),e(1),m(1),a(1),bits(30),i,n,ebit
      common /mpcomm/ junkc(600),junk1(120),junk2(120),junk3(120),junk4(
     x120)
      character junkc
      integer junk1,junk2,junk3,junk4
      call mpcopy(b,a)
      n=e(1)
23201 if(.not.(n.gt.1))goto 23203
      ebit=e(n)
      i=0
      do23204 j=1,30 
      i=i+1
      bits(i)=mod(ebit,2)
      ebit=ebit/2
      if(n.eq.e(1).and.ebit.eq.0)then
      i=i-1
      goto 23205
      endif
23204 continue
23205 continue
23208 if(.not.(i.gt.0))goto 23210
      call mul(a,a,junkc)
      call div(junkc,m,junk4,a)
      if(bits(i).ne.0)then
      call mul(a,b,junkc)
      call div(junkc,m,junk4,a)
      endif
23209 i=i-1
      goto 23208
23210 continue
23202 n=n-1
      goto 23201
23203 continue
      return
      end
      integer function jacobi(pa,qa)
      integer pa,qa,p,q,r,itmp
      character*80 e
      e = 'bad argument for JACOBI'//char(0)
      jacobi=1
      if(qa.le.0.or.pa.lt.0.or.mod(qa,2).eq.0)then
      call putdec(6,pa,5)
      call putdec(6,qa,5)
      call putc(char(10))
      call error(e)
      endif
      if(qa.eq.1.or.pa.eq.1)then
      return
      endif
      p=pa
      q=qa
      if(q.le.p)then
      p=mod(p,q)
      endif
      if(p.eq.0)then
      jacobi=0
      p=1
      endif
23221 if(p.ne.1)then
      itmp=1+and(p,3)
      I23223=(itmp)
      goto 23223
23225 continue
      p=p/4
      goto 23224
23226 continue
      r=mod(q,p)
      if(r.eq.0)then
      p=1
      jacobi=0
      else
      q=p
      p=r
      endif
      goto 23224
23229 continue
      p=p/2
      r=mod(q,8)
      if(r.eq.3.or.r.eq.5)then
      jacobi=-jacobi
      endif
      goto 23224
23232 continue
      if(mod(q,4).eq.3)then
      jacobi=-jacobi
      endif
      r=mod(q,p)
      if(r.eq.0)then
      p=1
      jacobi=0
      else
      q=p
      p=r
      endif
      goto 23224
23223 continue
      I23223=I232230
      if (I23223.lt.1.or.I23223.gt.4)goto 23224
      goto (23225,23226,23229,23232),I23223
23224 continue
      goto 23221
      endif
23222 continue
      return
      end
